#HKA tests
#first get genomic data for chicken
chicken.mart <- useMart('ensembl',dataset='ggallus_gene_ensembl')
tst <- getSequence(mart=chicken.mart,id="ENSGALT00000000024",type="ensembl_transcript_id",seqType="gene_exon_intron")
tst2 <- reverse(readEmblSeq('/Paterson/Datafiles/grouse/embl_exon3/Contig148_exon.embl'))
library(Biostrings)

tmp.align <- pairwiseAlignment(tst2,tst$gene_exon_intron,type='global')
#not really working

#try just using data from trimmed seqs
HKA.trimmed <- data.frame(grouse.id = grouse.chick.MK2$grouse.id,
	inheritance=rep(1,nrow(grouse.chick.MK2)),
	sp1.length=grouse.chick.MK2$aligned.length,
	sp2.length=grouse.chick.MK2$aligned.length,
	comp.length=grouse.chick.MK2$aligned.length,
	sp1.noseq=rep(100,nrow(grouse.chick.MK2)),
	sp2.noseq=rep(1,nrow(grouse.chick.MK2)),
	sp1.noP=rowSums(grouse.chick.MK2[,c('Pn','Ps')]),
	sp2.noP=rep(0,nrow(grouse.chick.MK2)),
	divergence=rowSums(grouse.chick.MK2[,c('Dn','Ds')]),
	stringsAsFactors=FALSE
	)

HKA.trimmed <- HKA.trimmed[HKA.trimmed$sp1.length>0,]

#get length of sp1 sequence properly
for(i in 1:nrow(HKA.trimmed)){
	tmp <- try(getCDS(dir="/Paterson/Datafiles/grouse/embl_exon3",seq.name=HKA.trimmed$grouse.id[i],single.cds=TRUE))
	if(class(tmp)=="try-error") next
	HKA.trimmed$sp1.length[i] <- nchar(tmp[1])
	}
HKA.trimmed$brief.locus <- paste("loc",1:nrow(HKA.trimmed),sep="")
#write out to file
write.table(HKA.trimmed[,c(11,2:10)],file="/Paterson/Datafiles/grouse/test_files/hka_trimmed.txt",sep="     ",row.names=FALSE,col.names=FALSE,quote=FALSE)
#hka doesn't like this

HKA.trimmed2 <- HKA.trimmed[HKA.trimmed$sp1.noP>0,]
write.table(HKA.trimmed2,file="/Paterson/Datafiles/grouse/test_files/hka_trimmed2.txt",sep="     ",row.names=FALSE,col.names=FALSE,quote=FALSE)
#still doesn't like it

write.table(HKA.trimmed2[1:50,],file="/Paterson/Datafiles/grouse/test_files/hka_trimmed_test.txt",sep="     ",row.names=FALSE,col.names=FALSE,quote=FALSE)


#then use snps from intronic seq???